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PREFACE 


This  Memorandvim  deals  with  a  new  mathematical 
approach  to  the  problan  of  handling  a  certain  type  of 
computationally  complex  equation  that  the  authors  have 
encountered  In  their  work  on  mathematical  models  for 
drug  distribution  In  the  body.  Such  studies  have 
Important  Implications  for  cancer  chemotherapy  as  well 
as  for  other  physiological  processes. 
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SUMMARY 

In  their  work  on  the  construction  of  mathematical 
models  to  Mill'd  dlfl’  the  study  of  physiological  processes 
connected  with  cancer  chemotherapy,  the  authors  nave 
encountered  convolution  equations  that  may  present- 
severe  computational  difficulties  due  to  the  storage 
requirements  for  their  solution,  In  this  Memorandum, 
the  authors  present  a  new  approach  to  this  problem, 
using  the  technique  of -differential  approximation. 
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DIFFERENTIAL  APPROXIMATION  APPLIED  TO  THE  SOLUTION 
OP  CONVOLUTION  EQUATIONS 

1.  INTRODUCTION 

In  the  course  of  constructing  some  mathematical 
models  of  physiological  processes  connected  vflth  cancer 
chemotherapy  tlJ,  vre  have  encountered  functional 
equations  containing  convolution  terms.  Equations  of 
this  type  are  unpleasant  computationally  because  of  the 
storage,  and  thus  time,  requirements  fpr  solution.  In 
some  cases,  these  storage  requirements  could  excee-d  .  - 

present  capabilities  and  thus  seriously  impede 
numerical  solution. 

We  v?ish  to  present  a  new  approach  to  this  problem 
using  the  technique  of  differential  approximation.  To 
illustrate  the  method,  we  shall  consider  the  equaM-oh.  - 

(l.l)  u(t)  »  f(t)  +  J"  u(s)ds. 

0 

2.  POLYNOMIAL  APPROXIMATION  AND  EXTENSIONS 

A  classical  problem,  vrfiich  owes  its  inception  to  a 
control  process  associated  with  the  Watt  steam  engine 
(see  [2]),  is  that  of  obtaining  a  polyriomial  which 
deviates  the  least  from  a  given  function,  where  the 
deviation  is  measured  by  an  assigned  nom. 

If  we  recognize  that  a  polynomial  Pj^(t)  *»  Sq  +  a^t 
+  •  *  •  +  a^t*^  is  a  solution  of  the  linear  differential 
equation 
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(2.1) 


d?^ 


0, 


then  we  see  Immediately  that  this  problem  is  a 
particular  case  of  the  more  general  problem  of  fincLLng 
an  equation  - 


d<"+l)u  .  .  d"u 


dt* 


+  a^(t)u 


0, 


^^Qse  solution  approximates  to  the  given  function  in  an 
optimal  fashion.  In  the  partieular  case  mier-e  the 
aji^(t)  are  constants,  this  Is  equivalent  to  asking  for 
approxlm.ation  by  an  exponential  polynomial 

N  X.  t 

(2.3)  p„(t)  -  S  qj^(t)e  ^  .  ■ 

“  lc=o  ^ 

This  problem  Is  in  turn  a  special  case  of  the 
general  problem  of  approximating  to  a  given  function 
f(t)  by  means  of  the  solution  of  a  nonlinear 
-  differential  equation 


(2.4) 


d^~~  ^  -  U 


g(u,uv,...,u^^^). 


This  Is  a  meaningful  approximation  problem  for 
arbitrary  n,  since  an  arbitrary  analytic  fmctlon  will 
not  In  general  satisfy  a  nonlinear  differential  equation 
of  any  finite  order;  e.g.,  r(t). 


';:sa:::^m)es!gsmsids^ 
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This  problem  arises  In  the  study  of  design  and 
control  and  has  an  important  role  In  the  study  of 
adaptive  processes  (see  isl#  IW).  We  Will  discuss 
these  matters  elsewhere. 

3.  LINEAR  DIFFERBStTIAL  APPROXIMATION 

We  wish  to  consider  the  problem  of  approximating  to 
a  given  fianctlon  f(t)  by  means  of  an  exponential 
polynomial  of  the  type  appearing  in  (2.3)>  for  reasons 
we  shall  describe  below.  Since  a  direct  approach  to 
this  problem  possesses  well^hnom  pitfalls  (see  Lanczos 
1 5]  ) ,  we  shall  pursue  a  different  path. 

First  of  all,  we  shall  suppose  that  the  given 
function  f (t)  satisfies  an  ordinary  differential 
equation  " 

(3‘1)  f^"'^  =  h(f,f',..  .,f^”^^\t), 

since  this  Is  quite  often  the  case  in  applications. 

Secondly,  we  shall  determine  an  approximating 
linear  differential  equation  with  constant  coefficients 

(3.2)  f^^^  +  aj^f^^^^  +  ajjf  =  G, 

by  asking  that  the  coefficients  a^^  be  chosen  so  as  to 
minimize  the  .functional 

(3.3)  f  (f^^^  +  ajf^^^^  +  ...  +  aj^)^dt, 

0 

where  f  is  determined  by  (3.1). 


M-aar'rt.'riaitii'tOTiltiif 


The  approximation  to  f  will  then  be  the  solution 
of  thb  linear  differential  equation 

(3.^)  u^^^  +  aj^u^^^^  +  ...  +  ajjU  »  0, 

with  initial  conditions  vdilch  will  be  determined  in  a 
fashion  discussed  below.  - 


4 .  COMPUTATIONAL  ASPECTS 

The  minimization  of  the  expression  in  (3.3)  leads 
to  the  system  of  N  simultaneous  linear  equations- 


(4.1) 


N  /  pT 


dt^ 


0,  4  »  o:,l,2,.../N 


For  moderate  values  of  N,  l.e.,  N  ^  20,  the  computa¬ 
tional  solution  provides  no  difficulty  once  we  have 
evaluated  the  integrals  appearing  as  coefficients*  W© 
covild,  if  we  so  desired,  integrate  by  parts  and  reduce  the 
evaluation  of  these  integrals  to  the  evaluation  of  the 

Integrals  (f^^^)^dt.  For  moderate  Size  N,  however, 
0 

it  is  more  convenient  to  proceed  directly  as  follows. 

0, 1,  ...  ,  N  —  1 , 

— u 

defined  by  the  equations 


Introduce  new  variables  u^^j,  i,J 


(4.2) 


du, 


^  «  f (N-*i)j.(j)^  ’“ij(O)  "  0> 


and  solve  these  equations  simultaneously  with  the 
original  equation  for  f,  namely. 
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{4.3)  fW  -  g(i-,f',...*f<*^^),t). 

5.  SOLUTION  OP  APPROXIMATE  LINEAR  EQUATION 

Having  determined  the  coefficients  by  means  of 
the  foregoing  procedures,  we  now  wish  to  determine  the 
function  u(t)  as  a  solution  of  {3«^)»  A  first  ^ 
approach  is  to  use  the  initial  values 

(5.1)  u^^^(O)  »  f^^^(O),  1  -  0,1,2, »..,N  -  1, 

and  indeed  this  is  what  we  do  below  with  some  success. 

In  general,  however,  we  would  proceed  in  the 
following  fashion.  Let  . . .  ,Ujj  be  the  N 

princlpail  solutions  of  (3.^),  the  solutions  deteimined 
by  the  condition  that  the  matrix  whose  coltmns  are 
(ui(0)#^2(o)#.  •  etc.,  is  the  identity  matrix. 

Every  solution  can  then  be  written  in  the  form 

(5.2)  u(t)  -  S-ovU^-r 

1=1  t  ^ 


vrtiere  the  c^^  are  scalars.  Let  us  choose  these  c^^ 
so  as  to  minimize  the  expression 


(5.3) 


r  T  /  N  n2 
0 


The  equations  for  the  c^  are 


pT 

(5.4)  J  fu^dt  - 
0 


N  P  T 

2  Cj  J  u^Ujdt  »■  0,  1 

3^1  0 


To  determine  the  various  integrals,  we  introduce  the 
variables  v^^  and  by  means  of  the  relations 


(5.5) 


dv 


^  »  fUj^,  v^(0)  **  0, 


dw. 


"3^ 


U 


UiUj,  w^j(O)  »  0, 


oin  the  equations  for  the  Uj^  (Eq.  (3.^)  with 


appropriate  boundary  conditions)  and  the  equation  for 
f,  and  integrate. 

6 .  AN  EXAMPLE 

Let  us  now  discuss  the  equation  of  renewal  type 
given  in  (l.l).  Talcing  0  t  ^  1,  we  obtain  as  a 
third-order  differential  approximation  to  the  fmctlon 
e  the  selutien  -sf--the  eqaation 

(6.1)  u^3^  +  2.740299u^^^  +  7.9511452u^^^ 

+  5.7636455U  =  0. 

-t^ 

Using  the  initial  values  obtained  from  e  ,  namely 


(6.2)  u{0)  =  1,  u'(0)  »  0,  u"(0)  «  -  2, 


we  found  such  excellent  agreement  with  the  values  of 
e  over  0  ^  t  ^  1  that  there  was  no  need  to  follow 
the  procedure  of  Sec.  5. 


I 


f 


•ii- 


y 


-7t- 

Consider  the  expression 

(6.3)  w(t)  •  J  k(t  —  s)u(s)ds. 

0 

Differentiating  repeatedly,  and  adding  with  the 
coefficients  obtained  above,  we  have 

(6.4)  +  2.740299w^^^  -t  7^9511^+52w^^^  +  5.7636455w 
-  k(0)u"(t)  +  k'(0)u»(t)  +  k"(0)u(t) 

+  2.740299[lc(0)u»(t)  +  k‘(0)u(t)] 

+  7.9511452  k(0)u{t) 

P  t  -  ' 

+  J  u(s)[fc»»'{t  -  s)  +  2. 7402991c" (t  ~  s) 

0 

+  7.95114521c‘ (t  -  s)  +  5.7636455k{t  ~  s)]  ds. 

Taking  k(t)  =  e  and  assuming  that  the  term  under 

the  Integral  sign  is  negligible,  we  obtain  a  third  order 

linear  differential  equation  for  w  =  u  —  f. 

Let  us  take  f  »  1  —  /^  e  ds  so  that  the  equation 

0 

(6.5)  u(t)  ’=>  1  -  e”®  dB  +  J  u(s)ds 

0  0 

has  the  solution  u(t)  =1. 

The  function  f(t)  as  given  above  satisfies  the 
linear  differential  equation 

(6.6)  f^^^  +  2tf^^^  +  2f^^^  =  0, 

with  f(0)  -1,  f'(0)  »  -  1,  f"(0)  «  0. 


-j 


Solving  (6., 6)  together  with  the  approximate  linear 
equation  for  u  obtained  from  (6.4)  j  we  obtain  the 
following  values  for  u(t): 


t 

u(t) 

u'(t) 

u"(t) 

0.1 

0.999999 

— 

0.2 

0.999999 

-0.146x10“3 

-0.148x10"^ 

0.3 

0.999969 

— 

0.4 

0.999937 

— 

0.5 

0.999909 

— 

— 

0.6 

0.999898 

— 

0.7 

0.999909 

0.229x10“^ 

0.174x10“3 

0.8 

0.999938 

0.330x10“^ 

0.167x10“"^ 

0,9 

0.999970- 

0.272x10“^ 

-0.135x10"^ 

1.0 

0.999989 

0.919x10“^ 

-0.189x10"^ 

As  we  can  see,  the  agreement  with  the  desired  value 
Ax(t)  =  1,  is  excellent^  ^-'-^-  . .  '  - -  = 

7.  DISCUSSION 

Consider  a  system  of  renewal-type  equations,  given, 
say, in  matrix  form  : 

(7.1)  X(t)  »  P(t)  K(t  -  s)X(s)ds. 

0 

Equations  of  this  type  arise  naturally  In  the  study  of 
multidlinensloneil  branching  processes,*  see  [5],  [6]. 

If  X(t)  Is  a  5  X  5  matrix,  we  are  required  to 
store  25  functions  (i.e.,  the  elements 
i,j  -  1,2,. ..,5)  if  we  proceed  in  th6-usual  fashion. 


0 


If  high  order  accuracy  were  required — say,  intervals  of 
lO""^  over  Q  t  ^  5~we  would  find  that  rapid-access 
storage  capacity  would  be  exceeded. 

On  the  other  hand,  if  we  use  the  foregoing  tech¬ 
nique,  differential  approximation  of  order  5  wotlld  lead 
to  the  task  of  solving  about  250  simultaneous 
differential  equations  plus  those  required  to  determine 
F(t)  .  This  is  a  simple  matter  for  a  modem  computer. 
Purthemore,  it  is  clear  that  we  could  use  an  approxi¬ 
mation  of  order  10  without  coming  close  to  the  storage 
capacity. 


BEPERENCES 


i  1.  Bellman,  R.,  J.  Jacquez,  R.  Kalaba,  and  B.  Kotkln, 

A  Mathematical  Model  of  Drug  Distribution  In  the 
Body;  Implications  for  Cancer  dhemobherapy,  'ihe 
Rand  corporation,  Rrt^346>-M1R,  February  1963. 

Jj 

}:  2.  Bellman,  R.,  and  R.  Kalaba,  Mathematical  Trends  in 

Control  Tiieo^,  Dover  Publications,  Mew  'X'ork,  to 
;  appear  in  19^5- 

3«  Bellman,  R.,  H.  Kaglwada,  and  R.  Kalaba,  A  Computa¬ 
tional  Procedure  for  Optimal  System  Design  and 
mTfzatTon,  The  RANIi  ^"^oratlon,'  RM-317^PR, 

June  1562";  also  published  in  Proc.  Nat.  Acad.  Sci. 
USA,  Vol.  43,  1962,  pp.  1524-IB55: 

>  4.  Bellman,  R. ,  "Mathanatical  Model-making  as  an 

Adaptive  Control  Process,"  Mathematical  Optimization 
Techniques.  University  of  California  Press, 

Berkeley,  1963,  PP.  333-339. 

5.  Harris,  T.  E.,  Branching  Processes,  Ergebnlsse  der 
Math.,  Springer,  Berlin,  to  appear  in  1963. 

>  6.  Beilman,  R. ,  and  K.  L.  Cooke,  Dirf exynt lal^^f f erenc  e 

f  Equations,  Academic  Press,  iJew  Yorlc,  1953. 


